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ABSTRACT 

Three-dimensional numerical magnetohydrodynamic (MHD) simulations are per- 
formed to investigate how a magnetically confined mountain on an accreting neutron 
star relaxes rcsistivcly. No evidence is found for non-ideal MHD instabilities on a 
short time-scale, such as the resistive ballooning mode or the tearing mode. Instead, 
the mountain relaxes gradually as matter is transported across magnetic surfaces on 
the diffusion time-scale, which evaluates to t\ ~ 10 5 — 10 8 yr (depending on the con- 
ductivity of the neutron star crust) for an accreted mass of M a = 1 .2 x 1O~ 4 M0. 
The magnetic dipolc moment simultaneously reemcrges as the screening currents dis- 
sipate over t\. For nonaxisymmctric mountains, ohmic dissipation tends to restore 
axisymmetry by magnetic reconnection at a filamentary neutral sheet in the equato- 
rial plane. Ideal-MHD oscillations on the Alfven time-scale, which can be excited by 
external influences, such as variations in the accretion torque, compress the magnetic 
field and hence decrease t\ by one order of magnitude relative to its standard value 
(as computed for the static configuration). The implications of long-lived mountains 
for gravitational wave emission from low-mass X-ray binaries are briefly explored. 
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1 INTRODUCTION 

Observations suggest that the magnetic dipole moment 
of accreting neutron stars i n X-ray binaries, fi, decreases 
with accreted mass, M a ( Taam fc van de Heuvell 1 19861 : 
Ivan den Heuvel fc Bitzarakil 19951). possibly through mag- 
netic screening or burial ( Bisnovat yi-Kogan fc Komberg 



1974 lRomanilfl990l : [Pay no 



Melatosl 12004 : lLovelace et al 



20051 ). During the burial process, the accreted plasma is 



channeled onto the magnetic poles of the neutron star, 
whence it spreads equa torwards, thereby distort ing the 
frozen- in magnetic flux (|Melatos fc Phinnevl 1200 ll ). Qua- 
sistatic sequences of ideal-magnetohydrodynamic (ideal- 
MHD) e quilibria describing how b urial proceeds were com- 
puted bv lPavne fc Melatosi (|2004 ) (hereafter PM04). These 
authors found that the magnetic field is compressed into an 
equatorial belt, which confines the accreted mountain at the 
poles. 

Surprisingly, the distorted equilibrium magnetic fields 
resulting from b urial are stable in ideal M HD. In an axisym- 
metric analysis, iPavne fc Melatosi (|2007l ) (hereafter PM07) 
found that the mountain, when perturbed, oscillates radi- 
ally and laterally in a superposition of global Alfven and 



compressional modes, but it remains intact. Of course, an 
axisymmetric analy s is neg lects important toroidal modes. 
IVigelius fc Melatosi (|2008rJ ) (hereafter VM08) found that 
the axisymmetric configuration is unstable to the undulat- 
ing submode of the three-dimensional Parker instability in 
spherical geometry. Again though, while the hydromagnetic 
structure reconfigures itself globally, the mountain remains 
confined to the magnetic poles once the instability saturates. 

PM04, PM07, and VM08 considered ideal-MHD equi- 
libria. However, magnetic burial creates steep magnetic 
gradients, which relax resistively. A conservative estimate 
of the relative importance of nonideal effects can be ar- 
rived at by assuming that the electrical resistivity in 
the outer crust is dominated by electron-phonon scatter- 
ing. Under this assumption, resistive relaxation arrests 
the growth of a mountain when the accreted mass ex- 
ceeds ~ 10~ 5 M W (iBrown fc Bildstenlll998l ; ICumming et al.l 
12004 iMelatos fc Payne! |2005| ).~ Resistive instabilities, like 
the global tearing mode or the local gravitational mode 
|Furth et alJll963h . grow faster than the simple resistive 
time-scale. Three-dimensional modes like the resistive bal- 
looning mode, which grows if the pressure gradient is parallel 
to the field line curvature, may rapidly destroy the confine- 
ment of the mountain. 
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During the early stages of accretion, the moun- 
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tain might be disrupted on the A lfven time s c ale b y 
the ideal-MHD balloonin g m ode dLitwin et all l200ll) . 
IVigelius fc Melatosl (|2008bh and IVieeliusI (|2008h show that 
equatorial magnetic stresses stabilize the configuration and 
generally prevent disruption in the high-M a regime. These 
authors then continue to solve the initial value problem by 
injecting plasma into an initially homogenuous background 
(with M a = 0) threaded by a dipolar field. They find no 
evidence for a growing instability in the low-M a regime. In 
this article, we investigate further how resistive relaxation 
competes with accretion at different accretion rates. 

The main aim of this article is to test if resistive in- 
stabilities disrupt the mountain on time-scales comparable 
to the accretion time-scale. The article is divided into six 
sections. Section [2] introduces the numerical setup used in 
our simulations, section [3] describes the dynamics of the re- 
sistive relaxation, and section [4] characterizes the magnetic 
field structure. In section [5] we evaluate the resistive relax- 
ation time as a function of accretion parameters. In section 
[6] we study how rapidly the magnetic field reemerges after 
accretion stops. We discuss our results in section]?] focussing 
on the ramifications for gravitational wave emission from ac- 
creting neutron stars. 



2 NUMERICAL MODEL 

2.1 Grid and units 

The simulations in thi s paper employ t he parallel, ideal- 
MHD solver zeus-mp jHaves et alj 120061 ). extended to in- 
clude resistive effects, as described in appendix [X] All the 
simulations are carried out in a spherical polar coordinate 
system (r,8,(f>), where r is logarithmically stretched as de- 
scribed in PM07 and VM08. To handle the disparate radial 
and lateral length-scales, we set up a downscaled neutron 
star with M, = 1.01 x 10~ 5 M Q and R, = 2.7 x 10 3 cm, such 
that the curvature a — R*/ho = 50 is still large while the 
hydrostatic scale height ho = 53.8 cm (defined in PM04) is 
preserved. We justify this approach by noting that the small- 
Ma analytic solution depends on M» and R, only through 
the combination ho ( PM04). The downscaling transforma- 
tion was employe d in iPavne fc Melatosl (120071) and VM08 
and validated by IVigelius fc Melatosl (|2008bF r in the large- 
Ma regime. 

Throughout this paper we fix pio — G — c s — ho = 1, 
such that the base units (in cgs) for mass, magnetic field, 
time, and resistivity become Mo = ho(? s /G = 8.1 x 10 24 
g, Bo = [moc 4 /(G^)] 1/2 = 7.2 x 10 17 G, r = h /c s = 
5.4 x 10~ 7 s, and 770 = Tq 1 = 1.86 x 10 6 s _1 respectively. 
The characteristic mass (PM04) then evaluates to M c = 

6.2 x 10~ 15 A/ Q for the downscaled star. 

2.2 Initial and boundary conditions 

Our aim in this paper is to examine the influence of a finite 
conductivity on magnetic mountain equilibria in 2.5 and 

3 dimensions. Axisymmetric equilibria are imported from 
the Grad-Shafranov (GS) solver developed by PM04. Non- 
axisymmetric equilibria are imported from zeus-mp after 
the transient, three-dimensional Parker instability saturates 
(VM08). All our simulations are isothermal (XIS0=. true . ). 



Boundary conditions are enforced in zeus-mp by ghost 
cells framing the active grid. Our choice of a spherical polar 
grid requires periodic boundary conditions at the <f> bound- 
aries [ikb.niks(l)=4 and ikb.noks(l)=4]. The 6 = tt/2 
boundary is reflecting [ojb.nojs(l)=5], with v± = By = 0. 
The line 9 = is also reflecting [ijb.nijs(l)= -1], with 
tangential magnetic field (v± = B± — 0). Additionally, the 
toroidal component B$ reverses at 6 = 0, i.e. B<f,(—6) = 
Bj,{6). The outer boundary at r = R m is a zero-gradient 
boundary [oib.nois(l)= 2]. The magnetic field at r = R* 
is line-tied by fixing the plasma variables [iib.niis(l)= 3] 
at this boundary: B is dipolar and p is kept several orders 
of magnitudes higher than the active grid values in order to 
realise an impenetrable surface. 



2.3 Resistivity 

The electrical conductivity a is a key input into the models 
presented in this article. In the outer crust, all transport pro- 
cesses are dominated by electrons scatter ing off phonons and 
impur ities [for a recent review compare IChamel fc Haensel 
|2008l )] and a can be derived from the scattering frequen- 
cies in the relaxation time approximation. For temperatures 
below the Umklapp temperature (Tu ~ 10 7 K), electron- 
phonon scattering is suppre ssed and the conductivity m ust 
be attributed to impurities (|Cumming et al.ll200ll . 120041 ') . In 
rapid accretors (M > 10 _11 M o yr" 1 ), one finds T > 10 s 
K and phonon scattering dominates, provided the impu- 
rity concentration satisfies Q < 1. In accreting neutron 
stars, Q is set by the composition of the ashes produced in 
steady state nuclear burning at low densities. ISchatz et al.l 
(| 199Sh find a large variety of nuclei in the crust so the im- 
purity factor is high (Q ~ 100). They argue that impurity 
scattering therefore dominates, except for very rapid accre- 
tors ( Q ~ 1 for M > 30M E dd). On the other hand. Ijonesl 
(|2004h noted that, if the primordial crust is completely re- 
placed by heterogeneous accreted matter, a temperature- 
independent conductivity dominates electron scattering and 
one finds Q J> 1. Most autho r s (|Ko nar fc Bhatta charval 
1 19971 ; ICumming et~"al1l200ll , 12004 iPons fc Geppertll2007l ) as 
sume Q -C 1, as do we. 

Neglecting impurities, IPotekhin et al.l (|l999f ) compute 
the frequency of electron-ion scattering in liquid and solid 
Fe matter for a v ariety of temperatures and densities. 
IChamel fc Haensell (|2008l ) present a computation of a i^n- 
cluding impurity scatterin g) for an accreted crust model 
|Haensel fc Zdunikl Il990ah , finding 23 «C log 10 (cr/s _1 ) 
27.4 for T = 10 7 K and 10 9 s: p[gcm" 3 ] < 10 13 (note that 
this density r ange covers the whole outer crust including 
neutron drip). ICumming et al.l |2004h find similar values for 
an accreted crust, viz. <r p = 1.8 x 10 25 s~ 1 (p 7 { 6 /T|) for 
electron-phonon and <tq = 4.4 x 10 25 s~ 1 pY i 3 for impurity 
scattering (provided Q = 1). 

To model a realistic star, we choose the electrical resis- 
tivity to be 1.3 x 10 -27 s < 7? r < 1-3 x 10~ 24 s, covering the 
range quoted in the previous paragraph. Throughout this 
paper, we also run simulations with artificially high values 
off), in the range 1.3xl0 -27 ^ (r?/ls) ^ 9.2 x 10" 11 , in order 
to accelerate resistive processes and observe their evolution 
over a computationally practical time interval. 

For simplicity, we assume an isothermal equation of 
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Table 1. Simulation parameters, rj measures the resistivity in 
terms of the realistic value r] T = 1.3 X 10 -27 s and also deter- 
mines the Lundquist number Lu = tfj/ta, i.e. the ratio of the 
resistive time-scale Trj to the Alfven time-scale t&_. Models A-D 
are axisymmetric; models E-H are nonaxisymmetric. All models 
are for M„ = M c . 



10 20 30 40 50 



Model 


login fa/'fc) 


Lu 


axisymmetric 


A 


1 


5.7 x 10 15 


yes 


B 


14.9 


8.0 x 10° 


yes 


C 


15.9 


8.0 x lO" 1 


yes 


D 


16.9 


8.0 x 10~ 3 


yes 


E 


1 


2.99 x 10 14 


no 


F 


14.9 


4.25 x 10 3 


no 


G 


15.9 


4.25 x 10 2 


no 


H 


16.9 


4.25 x 10 1 


no 



0.5 



u 0.0 




state throughout this article. During the late stages of ac- 
cretion (M a > 10 -3 Mq), however, the magnetic mountain 
mass is comparable to the mass of the neutron star crust and 
the model mountain contains a wide range of densities and 
temperatures as a function of depth. Pycnonuclear reactions 
in the deep regions (p > 10 12 g cm -3 ) feed thermal energy 
into an adiabatic mountain. The assumption of isothermal- 
ity breaks down and a realis tic equation of state for non- 
catalyzed matter is required l|Haensel fc Zduniklll990ah . In 
particular, the accreted material is expected to s olidify at 
densities > 10 g cm (|Haensel fc Zdunik|[l990bh and will 
sink i nto the crust, which need s to be modelled as an elastic 
solid (|Ushomirskv et al.l lioooh . In a self-consistent model, 
the electrical conductivity will be computed as a function 
of p and T. Furthermore, a strong magnetic field (B S> 10 9 
G) breaks the symmetry of electron transport proc esses and 
causes an anisotropic conductivity (|Potekhinll 19991 ) . The ef- 
fect of a realistic equation of state is subject of current work 
and the results will be presented elsewhere. 

In light of the discussion above, it is not immediately 
obvious at what location in the crust a needs to be eval- 
uated. In the end, however, we note that there are other 
deficiencies in our model which outweigh the uncertainties 
in the conductivity (most notably, sinking). In the context 
of this article, we therefore treat a as a fiducial parameter. 
In particular, we will show how the resistive relaxation time 
scales with a in section [5] 



3 RESISTIVE INSTABILITIES 

In general, MHD systems with a finite conductivity ex- 
hibit a plethora of resistive instabilities actin g on time-scales 
much shorter t han the diffusion time-scale (|Lifschit3ll989l ; 
lBiskamdll993h . Our first task is to find out if such insta- 
bilities are present here and on what time-scales they act. 
Table[T]lists the simulations performed to this end. We track 
the evolution of the mass ellipticity e as a convenient way 
to parametrize the evolution of the global hydromagnetic 
structure (VM08). 



Figure 1. Evolution of mass ellipticity e for different Lundquist 
numbers (from top to bottom) Lu = 5.7 X 10 15 ,8, 0.8, 8 X 10~ 3 
(solid, dotted, dashed, dash-dotted) for the axisymmetric models 
A— D with M a = M c , The time is measured in units of the Alfven 
time (bottom axis). The top axis measures the time in units of 
the diffusion time for model C. Clearly, e decays on the diffusive 
time-scale. 



3.1 Axisymmetric dynamics 

Fig. [T] displays e as a function of time for models A-D in 
table [U In order to find out if a genuine instability grows on 
an e-folding time-scale t\ < td , we artificially increase r\ and 
hence the Lundquist number Lu = td/ta (models B-D). 
Here, ta = Lp 1 ' 2 /B and td = L 2 a denote the Alfven and 
the diffusion time-scales, respectively, a is the electrical con- 
ductivity and L = (|B|/| V 2 B|) 1//2 is a characteristic length- 
scale. Clearly, L, ta, and td are functions of position and 
time. We minimize L and ta over the integration volume, 
finding (for the axisymmetric models) L — 1.05 X 10 _3 /in 
and t a = 105to respectively. 

During the first oscillation cycle in model B (Lu = 8), 
e declines more steeply than in model A before taper- 
ing off. This behaviour becomes more distinct in model C 
(Lu — 0.8), where e decreases rapidly, then plateaus when 
the equatorward motion of the mountain stops and subse- 
quently reverses. This cycle of decline followed by plateauing 
repeats several times while e tends to zero overall. 

Particularly interesting from a physical point of view is 
the behaviour of model D, with td <S ta- As rj is large, the 
magnetic field is unable to contain the mountain at the mag- 
netic pole. Consequently, the plasma slips through the field 
and falls towards the magnetic equator, where it is reflected 
at the boundary; that is, the mountain meets its counterpart 
centred at the other pole. As a result, e oscillates around the 
abscissa. A realistic neutron star never enters the regime 
td -C ta, but the tendency of the mountain to slip and 
bounce affects the dynamics for all values of td/ta, as dis- 
cussed in section HOI 

Fig. [2] shows the density contours (dashed curves) and 
projected magnetic flux surfaces (solid curves) for a merid- 
ional slice of model C. Snapshots are taken at £/ta = 
0,0.428,1.09,1.76,4.76,5.80. At t/r A = 0.428 and 1.76, e 
is in decline, according to Fig. [1] (model C, dashed line). At 
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Figure 2. Meridional section of model C at </ta = 0, 0.428, 1.09, 1.76, 4.76, 5.80 (top left to bottom right). Shown arc density contours 
(dashed curves) with values log 10 (p/ p' ) = —13, —12, —11, —10.7, —10.5, —10.3 and magnetic flux surfaces in cross-section (solid curves). 
The plasma diffuses through the flux surfaces while the magnetic field relaxes radially. 



t/rA = 0, 1.09 and 4.76, e is in a plateau. The configuration 
settles down at t = 5.80ta. 

The oscillations in Fig. [1] and Fig. [2] are driven by the 
hydrostatic pressure gradient perpendicular to the magnetic 
flux surfaces. Their amplitude remains bounded. Pressure- 
driven instabilities, such as the interchange or ballooning 
mode, grow when the field line curvature has a component 
along the pressure gradient (i.e. k- Vp > 0, where n = b-Vb 
and b — B/B) , a con figuration termed unfavourable curva- 
ture (|Lifschitj[l989l ). The top left panel of Fig. [2] shows 
clearly that the pressure gradient (which is proportional to 
the density gradient) in the ideal-MHD equilibrium is op- 
posed to the curvature, preventing the onset of a pressure- 
driven instability. Line tying also contributes to stability 
(VM08). 



3.2 Nonaxisymmetric dynamics 

The stability of an MHD system changes considerably upon 
passing from two to three dimensions. It turns out that, in 
the ideal case, the additional degree of freedom accomodates 
toroidal Parker modes that rearrange the axisymmetric equi- 
librium into a slightly nonaxisymmetric state (VM08). The 
stability of this state when resistivity is switched on is the 
concern of this section. The relevant models are labelled E- 
H in tabled] 

Following section 13.11 we first examine the time evolu- 
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Figure 3. Evolution of mass ellipticity e for different Lundquist 
numbers (from top to bottom) Lu = 2.99 X 10 14 , 4.25 X 10 3 , 4.25 X 
10 ,42.5 (solid, dotted, dashed, dash-dotted) for the nonaxisym- 
metric models E-H with M a = M c . The time is measured in units 
of the Alfven time (bottom axis). The top axis measures the time 
in units of the diffusion time for model G. As in Fig. [T] t decays 
on the diffusive time-scale. 

tion of e for models E-H. The results are summarized in 
Fig. [3] Strictly speaking, the definition of e is only meaning- 
ful for an axisymmetric configuration. However, the three- 
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Figure 5. Mass cllipticity e (solid curve, left linear axis) and 
diffusion time-scale tq (dashed curve, right logarithmic axis) as 
functions of time, in units of the Alfven time, for model A. Both 
quantities are normalized to their initial values. During the first 
cycle, Trj drops by 91 per cent. 
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Figure 6. Temporal evolution of the total, magnetic, gravita- 
tional, kinetic, and acoustic energies W, W m , W g , Wk, and W a 
(top to bottom) for model C, all normalised to Wq = 2.4 X 10 36 
erg and corrected for mass loss through the outer border. Note 
the decrease in W m due to magnetic dissipation. 



dimensional equilibrium deviates from axisymmetry by less 
than 0.8 per cent (VM08), so e is a good proxy for the global 
hydromagnetic structure. We find that Model E is stable for 
t < 3.5ta- In models F and G, which have Lu < 4.25 x 10 3 
and ta = 124Yo, the mountain dissipates on the diffu- 
sive time-scale (e.g. td = 5.26to for model G). Model H 
(Lu — 42.5) exhibits the pressure-driven oscillations ob- 
served in model D (cf. Fig. [1]). 

The three-dimensional hydromagnetic structure of 
model G is captured in a series of snapshots in Fig. U 
Shown is the mountain (orange surface), delineated by the 
isosurface p = 1.03 x 10 9 g cm" 3 , along with the mag- 
netic field lines (blue and green curves), at the instants 
t/TA = 0,1.70,3.40,5.09,6.79,8.57. The initial configura- 
tion (top-left panel) is the outcome of the three-dimensional 
undulating submode of the Parker instability (VM08). The 
field lines curve towards the magnetic poles, while the or- 
ange isosurface spreads equatorwards by 32 per cent relative 
to its initial position. Soon after the resistivity is switched 
on (top-middle panel), the system behaves like model C: 
magnetic tension straightens the field lines radially, while 
the plasma slips laterally through the flux surfaces, allow- 
ing the magnetic mountain to escape its polar confinement 
and spread over the neutron star surface. However, the non- 
axisymmetric configuration is the saturation state of the 
transient Parker instability. Hence, unlike model C, the 
global hydromagnetic oscillations in model G have already 
died away. The instability time-scale is given by the diffu- 
sion time-scale, no t the tearing-mode time-scale (tdta) 1 ^ 2 
(|Furth et alJll963h . 



3.3 Oscillation enhanced diffusion 

As the axisymmetric mountain oscillates laterally, the field 
gradients steepen whenever the field compresses. This effect 
accelerates resistive relaxation. Fig. [S] plots td (right, loga- 
rithmic axis) and e (left, linear axis) as functions of time for 
model A. During the first cycle, td drops to nine per cent 



of its original value and diffusion proceeds proportionally 
faster. The effect of diffusion is two-fold, (i) The plasma slips 
through magnetic flux surfaces and moves towards the mag- 
netic equator. Eventually, as seen in the lower middle panel 
in Fig. [21 it covers the surface evenly and e decreases (Fig. 
[1]). (ii) Magnetic tension causes the field lines to straighten 
radially. Close to the magnetic equator, the hydrostatic pres- 
sure from the drained plasma also drives the magnetic field 
out wards. 

iMouschoviasI (1 19741 ) showed that an isothermal gravi- 
tating MHD system possesses a total energy W , which can 
be written as the sum of gravitational (W s ), kinetic (Wk), 
magnetic (Wm), and acoustic (Wa,) contributions, defined by 
Eqs. (10)-(13) in VM08. In ideal MHD, W is a conserved 
quantity. Adding resistivity allows the magnetic flux to dis- 
sipate, converting W m to W a via a source term (7 — l)?7|j| 2 
in the energy equation, where 7 is the adiabatic index. In an 
isothermal setup, this source term vanishes and the energy 
equation is trivially satisfied; heat is absorbed by a reservoir. 

The time dependence of the above four contributions to 
the energy integral for the axisymmetric model C are shown 
in Fig. [6] Following VM08, we correct for mass loss through 
the r = R m border by multiplying W g , Wk, and W a by 
M(t = 0)/M(t), where M(t) is the total mass in the simula- 
tion volume at time t. Clearly, some energy is converted to 
heat: W drops by 1.4 per cent during the interval t < 3ta, 
as the magnetic field dissipates. W g decreases because the 
accreted matter, which is initially confined at the magnetic 
pole, distributes itself evenly over the star's surface. W^ rises 
sharply when the whole system reconfigures and then slowly 
decreases due to numerical dissipation. Wa decreases along 
with |Vp|. 

Fig. 0shows the time dependence of the different energy 
contributions for the nonaxisymmetric model G. Similar to 
Fig. [6] W drops by ~ 2 per cent on the diffusion time- 
scale. The main losses occur in W m , which drops by one 
order of magnitude, and W g , which decreases by 1 per cent 
(from a high base). The kinetic energy slowly rises, as an 
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Figure 4. Density and magnetic field structure of model G at t/t\ = 0, 1.70, 3.40, 5.09, 6.79, 8.57 (from top left to bottom right). The 
mountain (orange surface) is defined by the isosurface p(r, 9, <ft) = 1.03 X 10 9 g cm -3 . In order to assist with visualization, all length-scales 
of the mountain and the field lines are magnified five-fold. The footpoints of the blue field lines start from the stellar surface, while green 
field lines start from the equatorial plane. 
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Figure 7. The evolution of total, magnetic, gravitational, kinetic, 
and acoustic energies W, Wm, Wg, Wk, and W a (top to bottom) 
for model G, all normalised to Wq = 2.1 X 10 36 erg, as a function 
of time (in units of the Alfven time) and corrected for mass loss 
through the outer border. 



overstable mode grows (see sect ion [6~2]l . Since zeus-mp does 
not explicitly include viscosity, Wk dissipates numerically 
(i.e. through the grid viscosity). 

On the other hand, discretizing the continuous MHD 
equations introduces numerical errors that dissipate mag- 
netic energy and can damp the growth of unstable modes. 
This numerical viscosity, u, can therefore artificially stabilize 
our configuration. A good measure for the relative impor- 
tance of v is the magnetic Prandtl number, Pr m = td /t v [ sc , 
where r v i sc = pl^/v with l v being a characteristic length 
scale for velocity gradients and v the viscosity. Since v 



owes its existence to the discretization of the MHD equa- 
tions, it depends on the grid size and the field gradients. 
In order to obtain an accurate estimate for Pr m , we com- 
pute the timescale, r v i ac , on which ideal- M HD oscillations of 
an ax isymmetric configuration die away (|Pavne fe Melatosl 
120061 ) and compare it to the diffusive timescale, Tb, finding 
Pr m = 1.03Z/U. Hence, the contribution of numerical vis- 
cosity is generally small (e.g. Pr m ~ 1CP 3 C 1 in model 
D). 



4 MAGNETIC FIELD STRUCTURE 

The global hydromagnetic evolution observed in section [3] 
occurs on the ohmic time-scale. This indicates that relax- 
ation is dictated by magnetic diffusion rather than resis- 
tive transient instabilities on short time-scales, such as the 
large-scale tearing mode or the localized gravitational mode 
l|Furth et al.ll 19631 ). Transient instabilites occur in the neigh- 
bourhood of current sheets, which dissolve into magnetic 
islands and dissipate. In this section, we examine the mag- 
netic geometry of the resistively relaxing mountain to check 
whether it is consistent with the above view that diffusion 
on large scales dominates the evolution. 



4.1 Neutral surfaces 

We begin by investigating the magnetic field structure of 
the axisymmetric model C (Fig. [2}. The initial equilibrium 
configuration is depicted in the top-left panel. Notice that 
there are no magnetic neutral points present. The mountain 
is held in place by the tension of the line-tied magnetic field. 

Does the co nfiguration contain current sheets? 
lHanasz et ail (|2002l ) showed that the undulating submode 
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Figure 8. Equatorial slices of the magnetic structure and current flows in model G at at = 0, 1.69, 3.39, 5.08, 6.77, 8.55 (from top 

left to bottom right). Shown are the projections of the magnetic field vectors onto the plane 8 = 1.5 rad (arrows) and the modulus of 
the current density \j\ (color coded). 



of the Parker instability in a Cartesian geometry creates 
current sheets in the plane perpendicular to the magnetic 
flux surfaces between regions with alternating polarity. Fig. 
[5] displays a time series of equatorial slices of the magnetic 
field from model G at t/r A = 0, 1.69, 3.39, 5.08, 6.77, 8.55. 
The projection of B/B onto the equatorial plane is in- 
dicated by arrows, while the current density |j| is color 
coded. The top left panel shows the initial configuration for 
our experiment, generated from an axisymmetric mountain 
after the undulating submode of the Parker instability 
saturates. While \j\ is greatest close to the stellar surface, 
where B is high, long radial current filam ents are also 
clearly present, albeit not as distinctly as in lHanasz et alj 
|2002l ). The filaments are neutral sheets. 



4.2 Reconnection 

Reconnection occurs at the current sheets in Fig. [8] quickly 
smoothing the toroidal gradients. Line tying at the stellar 
surface forces the field lines to adjust into a dipolar config- 
uration. A finite resistivity therefore acts to restore axisym- 
metry. In addition, the line-tying boundary condition acts 
as a source of magnetic flux, which is thence transported 

radially outward by diffusion. 

Where does reconnection occur? ISchindler" et al.1 (|l988l ) 
pointed out that a necessary and sufficient condition for 



■ 




Figure 9. Meridional slice of the local field-aligned electric field 
log 10 |£MJ|, normalised to BqEq = 1.51 X 10 _11 .Bo/Iq for model 
G. Reconnection mainly occurs near the stellar surface in the 
magnetic belt but also in the outer equatorial region (at r > 56/io 
and 0.7 <6 < 1.5). 
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global magnetic reconnection along some field line C is that 
the electric field has a component parallel to B, 

As E ■ B 0, (1) 

Jo 

where the integral is taken along C. (Equivalently, the helic- 
ity changes with time.) We plot a meridional slice of E ■ B 
at (j) = 2.3 rad in Fig. [9] Not surprisingly, E ■ B is high- 
est in the magnetic belt region, close to the star's surface. 
However, the undulating submode of the Parker instability 
also induces small toroidal currents (top left panel in Fig. 
[8|, so that E ■ B is high in the equatorial region too. We 
integrate E ■ B along two sample field lines with footpoints 
at (x ,e ,4> ) = (0,0.2,4.03) (field line ©) and (0,1.0,4.03) 
(field line ©) and find J c E ■ B = -8.6 x 10 _14 B (field 
line ©) and J c E ■ B = -3.2 x 10~ 16 B (field line ©) re- 
spectively. The topology of the magnetic field is discussed 
in section 14.31 where we show that field line CD undergoes 
reconnection while field line © does not. 



4.3 Topology 

In this subsection, we briefly discuss the change in magnetic 
topology brought about by reconnection. Fig. llOl displavs the 
magnetic field lines (solid curves) in three meridional slices 
= 0, 1.96, 3.93 (left, middle, right columns) for model G. 
One immediately notices that there is a Y-point located at 
(r,9) « (56/io, 1.4rad) in the top-left panel of the figure. The 
Y-point owes its existence to a boundary effect in the ideal- 
MHD simulation: during the onset of the Parker instability, 
the plasma is pushed out of the integration volume through 
the outer boundary. The subsequent backflow topologically 
separates the previously connected field lines. 

Associated with the Y-point is a current sheet at 9 ~ 1.4 
rad, which meanders like a band in the cj> direction. A cur- 
rent sheet naturally triggers reconnection. The bottom row 
of Fig. [TOl shows the same slices as the top row after 0.6td- 
Indeed, the field lines have reconnected: they are not topo- 
logically separated anymore, and the current sheet has van- 
ished. Alternatively, it is conceivable that the current sheet 
moves along with the plasma flow from its initial position at 
9 — 1.4 to the upper boundary at 9 — n/2. 

The concept of rational magnetic surfaces, where the 
field lines close upon themselve s, plays an imp ortant role in a 
local plasma stability analysis l|Lifschitj| 19891 ) . The bending 
of field lines as a result of a Lagrangian displacement £ is 
associated with an increase in potential energy, given by 
B • V£. In a tokamak geometry, it can be shown that this 
term vanishes on a rational surface, which is directly related 
to the pitch angle B^/B p : the safety factor is defined as 
q = AB^/ABp, and a rational surface is one where q is a 
rational number. In Fig. 1111 we plot the pitch angle as a 
function of the arc-length coordinate r\ (right panels) for 
four different field lines (labelled ©-© in the left panel) in 
model G. Close to the pole (lines ©-®), the pitch angle 
stays below ~ 3 per cent. For line ©, it increases towards 
the equator, ultimately reaching w 20 per cent. The zero 
crossing for © indicates that B$ changes sign, a relic of the 
undulating submode of the Parker instability which gives 
birth to this state. Diffusion does not eliminate the toroidal 
component completely. 



We attempted for completeness to characterize the 
magnetic topology near the neutral sur face using scale in- 
variants of the strain tensor dB j /dxj (|Chong et al.l Il990l : 
IParnell et al]|l996l : IPeralta et all 120081 ). but this approach 
yields ambiguous results in this instance. 



5 RELAXATION TIME 

We are now in a position to compute how long it takes for a 
magnetically confined mountain to relax resistively, given 77 
and M a - Ultimately, as t — > 00, the mountain spreads itself 
uniformly over the stellar surface [i.e. p — p(r)], threaded by 
a dipole field (i.e. j = everywhere). However, this process 
does not approach completion for realistic rj over the lifetime 
of an accreting neutron star. 

Let us define the ohmic relaxation time to be the time 
that elapses before the mountain relaxes to e _1 its initial 
cllipticity. Fig.[T]presents e(t) for an axisymmetric mountain 
with M a = M c as a function of the conductivity a = r\~ x . 
Reading off n from e(-n) = e _1 e(0), and fitting the trend by 
linear least squares, we obtain 

— = 1.7 x W~ 3 a, (2) 

TO 

where a is measured in units of 00 = ^o" 1 = 1-86 X 10 6 s _1 . 
For the upscaled star with a realistic i] T we find t\ = 6.3 x 10 6 
yr, which is comparable to the fiducial accretion time-scale 
race = 10 6 - 10 7 yr. 

Fig. [3] presents e(t) for a nonaxisymmetric mountain 
with M a = M c . Applying the same procedure from the pre- 
vious paragraph to Fig. [3] we find 

— = 0.02a. (3) 

n> 

For an upscaled neutron star with realistic rj T , we find n = 
7.6 x 10 7 yr, comparable to the fiducial accretion time-scale. 
Astrophysically, this is the key result of this paper: magnetic 
mountains in three dimensions relax resistively over ~ 10 5 — 
10 8 yr, (depending on the particular value of a; see section 
I2.3|l . not over shorter time-scales like ta and (tatd) 1//2 . Note 
that ro is a local quantity for a stationary mountain; n is a 
better measure of the global diffusion time. 

We compare r\ to the growth time of the resistive Parker 
instability, whose dispersion relation is calculated in ap- 
pendix [B] The growth time is shortest for short- wavelength 
modes and is therefore set by the grid scale (k ~ 6.6/i ~ 1 ) in 
our units. Also, the ratio of magnetic pressure to gas pres- 
sure, a, is maximal in the magnetic belt region, where the 
magnetic pressure balances the gas pressure, viz. a ~ 1. 
The growth rate is independent of k and is given by T — 
[1/2(1 - 2q) 2 ][q/(1 + a)](g 2 /k 2 u 2 )(i/-m), where r D is the 
time required to diffuse over one scale height. Applying equa- 
tion (|B17[) to the axisymmetric model C, we find the Parker 
growth-rate to be F = 4.07 x 10 6 r _1 . Fig. Q] shows clearly 
that t\ 3> r _1 , further supporting our conclusion that the 
resistive relaxation occurs on the diffusion time-scale and 
does not involve MHD instabilities. The same conclusion 
applies for the nonaxisymmetric model G, with the same 
growth rate as for the axisymmetric model. 

Another way to present the results on n is to ask how 
e varies with the accretion rate M. The simulations under- 
lying Fig. [12] differ from the others in this paper in one 
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Figure 10. Density contours (dashed) and magnetic field lines (solid) for meridional slices at rj> = 0, 1.96, 3.93 rad (left to right columns) 
and t/ro = 0, 201 (top and bottom row) for model G. The direction of the magnetic field is indicated by arrows. At the Y-point in the 
top-right corner of each panel, reconnection occurs. 




Figure 11. Magnetic pitch angle B^/Bp (right panel) as a function of arc length rj along four magnetic field lines (D-@ for model G, for 
a snapshot taken at t = 0.6i"rj. Field lines (blue curves) are identified in the left panel. The mountain is defined by the orange isosurfacc 
p(r, 9, 4>) = 1.04 X 10 9 g cm 3 . Red indicates the neutron star surface r = R*. In order to help visualize the structure, all length scales of 
the mountain and the field lines are magnified five-fold. 
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Figure 12. Mass ellipticity e versus accretion rate M (in units of 
10~ 4 Mp) yr -1 ) for an axisy mmetric grown mountain [see text and 
I Vigelius fc Melatosi ll2008bt) 1. The solid (dashed) curve represents 
a mountain aged t = 10 6 yr (t = 10 s yr). Although e is generally 
higher for the older mountain, since more mass has been accreted, 
it almost touches the curve for the younger mountain, because the 
resistive instability acts to reduce e. 

important respect: the mountain is grown from scratch over 
time (starting from M a = 0), with mass injected at the poles 
of an initially dipolar magnetic field, at a rate M and with 
77 7^ throughout the experiment. In other words, resis- 
tive relaxation competes simultaneously with accretion. By 
contrast, in Figs. [H 411I a Grad-Shafranov equilibrium is im- 
ported into zeus-mp, 77 is switched on at t = 0, and the 
mountain subsequently relaxes. Growing the mountain con- 
fers several advantages: it reflects the astrophysical process 
of burial more faithfully and enables us to reach M a = 10M C , 
cf. M a < 1AM C with the Grad-Shafranov method. The dis- 
advantage is that, at present, we cannot study how the 
mountain relaxes after accretion stops, because ZEUS-MP 
fails when the injection "nozzles" are turned off (suddenly 
or with taper), due to a numerical instability (the grown 
mountain contains nonzero flows). We are therefore unable 
to compare the two numerical experiments exactly, although 
they are in close qualitative agreement. A detailed explana- 
tion of t he injection algorithm and v erification tests can be 
found in I Vigelius fc Melatosi (|2008bh . 

A crucial question is whether the Grad-Shafranov equi- 
libria can be uniquely attained as accretion onto the mag- 
netic poles occurs, in particular, when 77 7^ 0. The ex- 
periments co nducted by growing the m o untain a b initi o 
[Fig. HI and IVigelius fc Melatosi (|2008br ): IVigeliusI <|2008h ] 
mimic time-dependent accretion more faithfully. The in- 
falling plasma continuously deforms an initially dipolar field 
and every snapshot represents the equilibrium configuration 
for a particular M a . These equilibria are in good agree- 
ment with previous results obtained analytically or nu- 
merically with the Grad-S hafranov code [cf. Fig. 4.11 in 
IVigelius fc Melatosi (|2008bl )]. In particular, we find no ev- 
idence for (ideal or resistive) instabilities occuring in the 
low-M a regime. 

On the other hand, the unavoidably finite size of the 
simulation box leads to a subtle uniqueness problem. The 
material that is added to the pole pushes the field lines to- 



wards the equator. Because of the boundary conditions we 
use, these field lines jump discontinuously when they touch 
the bottom right-hand corner of the box from d r B = 
when penetrating the boundary r — R m to B r = when 
penetrating the boundary 6 — n/2 (compare the top-right 
corner of the top-right and bottom-left panels in Fig. 1 10 p . In 
effect, this is a "reconnection-type" event which changes the 
topology of the field lines, their connectivity to the "out- 
side world", and therefore the effective functional form of 
dM/dtp (which we assume to be constant throughout the 
run). In practice, it is likely the effect is very small, the 
evidence being (i) the small mass outflow (< 1 per cent 
of the total mass) through r = R m during a typical run, 
and (ii) the very similar equilibria obtained from solving 
the Gr ad-Shafranov equation and growing the mountain ab 
initio l| Vigelius fc Melatosi |2008bh . In principle, though, it 
can lead to different final states if the mountain is grown 
with and without resistivity turned orQ. 

Fig. [T2] displays the ellipticity as a function of M for 
a young (t = 10 6 yr, solid curve) and an old (t = 10 s yr, 
dashed curve) object. To perform the simulation over a prac- 
tical length of time, we artificially increase 77 to 10 13,85 77 r 
(solid curve) and 10 15 85 77 r (dashed curve). We then use the 
scaling n oc 77" 1 derived from Figs. Q] and [3] to relate the 
results to astrophysical time-scales. We point out that each 
curve in Fig. [12] basically displays e(t) and we can relabel 
the abscissa using M a = Mt. 

There are two opposing effects in the figure. First, the 
older object has generally higher e for a given M, simply 
because M a is higher. Second, resistive relaxation has more 
time to reduce e in the older object, so the two curves almost 
touch at M » 5 x 1O -8 M0 yr -1 . The injection algorithm 
induces global hydromagnetic perturbations; these numeri- 
cal artifacts are visible as oscillations at the high-M end of 
either curve. 

How does our relaxation time compare to previ o us es - 
timates? In the small-M a regime, IMelatos fc Payne! (|2005l ) 
found analytically that resistive relaxation stalls mountain 
growth at e ~ 10 -5 (assuming electron-phonon scattering 
with a crustal temperature of T = 10 8 K). Our results sug- 
gest that a mountain with M a — M c = 1.2 x 1CP 4 M@ relaxes 
resistively over ~ 10 5 — 10 s yr. Furthermore, when accre- 
tion and relaxation proceed together, we find again that e 
saturates at ~ 1 ~ 5 , e ven for M a > M c , in accord with 
IMelatos fc Pavnel (|2005h . 

Similar estimates were given by iBrown fc Bildstenl 

(|l998l ) who evaluated the diffusion time in the crust. Taking 
into account electron-phonon and electron-impurity scatter- 
ing, they found that phonon scattering dominates impu- 
rity scattering (provided Q < 1) and td ~ 10 4 yr when 
the star accretes at the Eddington rate. However, these au- 
thors considered only spherically symmetr ic accretion and 
disreg arded the global magnetic structure. ICumming et al.l 
|200J) found td ~ 10 8 yr for a crustal temperature of 
T — 10 6 K, in accord with our results. 

We conclude this subsection with a brief discussion of 
how the relaxation time changes with M a . Fig.Q3]compares e 
for two axisymmetric models with different accreted masses, 
M a = 0.6M C (solid curve) and M a = 1AM C (dashed curve), 



1 Sterl Phinney, private communication 



Resistive relaxation of a magnetically confined mountain on an accreting neutron star 11 




o.o r 

2 3 4 5 6 

t/T B 



Figure 13. Evolution of mass ellipticity e as a function of time in 
units of the respective diffusion times To for axisymmetric models 
with M a = 0.6M C (solid curve) and M a = 1.4M C (dashed curve). 
The Lundquist number for both models is Lu = 1. The diffusion 
time is t d = 123to (t d = 66.4r ) for M a = 0.6M C (M a = 1.4M C ). 
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Figure 14. Normalized magnetic dipole moment dio/Rf^ versus 
time (in units of the diffusion time) for model C. The resisitivc 
instability relaxes the magnetic field radially, such that dio peaks 
at t = 7ttj. Eventually, the screening currents in the mountain 
dissipate and dig tends to its initial value. 



but the same Lundquist number Lu — 1. Note that e is 
displayed as a function of time in units of the respective 
diffusion time, r D = 123r (td = 66.4r ) for M a = 0.6M C 
{M a = 1.4M C ). The magnetic field of the M a = 1.4M C model 
is more distorted and, consequently, td is shorter. Both mod- 
els exhibit resistive relaxation on the diffusion time-scale. 



6 REEMERGENCE OF THE BURIED 
MAGNETIC FIELD 

6.1 Magnetic dipole moment 

An important diagnostic of the global magnetic struc- 
ture is its magnetic dipole moment. This integrated value 
has the advantage that it is o bservationally accessible 
( van den Heuvel fc Bitzaraklll995T l. Indeed, the observed re- 
duction of the magnetic dipole moment by accretion is a key 
motivation of the magnetic mountain concept (PM04). 

Following VM08, we define the magnetic multipole mo- 
ment tensor as 

d ij (r)=r i+1 J dQY*r-B, (4) 

where Yij denotes the spherical harmonics and r is the po- 
sition vector. Henceforth, we evaluate dij at the simulation 
boundary r = R m and drop r. 

The evolution of dio, plotted in Fig. 1141 illustrates the 
effect of ohmic diffusion on the magnetic structure. Initially, 
dio is buried by the distorted magnetic field. The resistive in- 
stability then allows B to straighten radially and reduce the 
field line curvature (cf. bottom-middle panel in Fig. , as 
described in section [3~T1 Ultimately, the line- tying condition 
of the inner boundary forces dio to approach the underlying 
dipole moment of the star before accretion, as the screening 
currents in the mountain dissipate. In this sense, one can 
say that the buried magnetic field reemerges. 

The physical mechanism behind reemergence is illumi- 
nated by examining the radial dependence of the dipole mo- 
ment, snapshots of which are plotted in Fig. 1151 plotted at 
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Figure 15. Magnetic dipole moment dio(r) [normalised to 
dio{r = 0)] versus radius (in units of ho) for model C at times 
t/rpj = 0,333,665 (solid, dotted, dashed curves). 



t/ro = 0, 333, 665 (solid, dotted, and dashed curves, respec- 
tively) . Initially, dio is screened within a thin layer near the 
surface. As the screening currents dissipate resistively, mag- 
netic flux is transported radially outward, thereby increasing 
the dipole moment measured by an outside observer. 

The nonvanishing components of the magnetic dipole 
and quadrupole tensors are displayed in Fig. [16] for an non- 
axisymmetric mountain (model G) . As for the axisymmetric 
case (Fig. I14|l , dio (bottom panel) increases over the diffu- 
sion time-scale as the magnetic field relaxes, tending to the 
underlying value at r = R*. Overall, dio varies by less than 
10 per cent over the simulation. The magnetic quadrupole 
moment d2i approaches zero as diffusion restores the dipolar 
field. 
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Figure 16. Magnetic dipole moment dio/R^ and magnetic 
quadrupolc moment d.2\IR% l for model G, normalised to the ini- 
tial value of dio/-R^j = 7.3 X 10 11 G, as a function of time (in 
units of the diffusion time scale td). All other components of the 
tensor d{j vanish due to symmetry. 




-0.4 t , , , J 

50 100 150 200 

t/T 



Figure 17. Components of the mass quadrupole moment tensor 
normalised to the maximum of Q33, 7.2x 10 24 g cm 2 , as a function 
of time, in units of the diffusion time-scale, for model G. 

6.2 TViaxiality 

The distorted magnetic field structure in Fig. [2] is accompa- 
nied by deformation of the mass distribution. This matters 
when considering accreting neutron stars as gravitational 
wave sources. Fig. [T7] plots the components of the Cartesian 
mass quadrupole moment, defined as 

Qij = J d s x (3x'iXj - r' 2 5 lJ )p{x'), (5) 

versus time for the nonaxisymmetric model G. The diagonal 
elements of Qij measure the axisymmetric distortion and are 
directly related to the ellipticity by e a Q22 oc Q33. As the 
mountain relaxes resistively, Q22 and Q33 decrease, asymp- 
toting at ~ 20 per cent of the initial value at t ~ 150td- 
This is normal: plasma diffuses across the flux surfaces and 
spreads evenly over the stellar surface. However, since mag- 
netic diffusion tends to smooth out gradients in B, thereby 



reducing the actual diffusion time-scale, we reach an inter- 
mediate, metastable state. In this state, the field lines are 
almost radial but the mountain has not yet diffused to cover 
the surface evenly. We compute the diffusion time-scale of 
the metastable state to be r{y = 27.6ro, five times higher 
than td = 5.26to of the initial state. Eventually, the remain- 
ing plasma diffuses over the time-scale t d and Qij tends to 
zero. 

The offdiagonal elements of Qij (top panels of Fig. I17|l 
measure the deviation from axisymmetry. They decrease on 
the time-scale n and then oscillate around the abscissa. Non- 
axisymmetric oscillations, observed previously in ideal-MHD 
calculations (VM08), are excited here when the mountain re- 
configures: small numerical inaccuracies perturb the steady- 
state equilibrium and the mountain readjusts on the Alfven 
time-scale. The period is ~ 10td for model G. The am- 
plitude initially grows then decays. The existence of such 
overstable modes is peculiar to a dissipative MHD system. 
The linear force operator is no longer self-adjoint and its 
eigenvalues generally have both a real and an imaginary 
part. The tendency of Q12 to decrease tallies with the obser- 
vation that resistivity restores axisymmetry by smoothing 
toroidal gradients, as postulated in section [47T1 It is impor- 
tant to note that the observed oscillations cannot arise if 
there is a realistic separation between the Alfven and diffu- 
sion timescale. For completeness, we note that the number 
Q11 — Q22I = IQ33 + 2Q22I is another measure for the de- 
parture from axisymmetry. However, it is obvious from Fig. 
[17] that the magnitude of this number is small compared to 
the magnitude of the diagonal elements. 



7 DISCUSSION 

The formation of magnetically confined mountains at the 
poles of an accreting neutron star is one explanation of the 
observed reduction of the magnetic dipole moment with M a 
in neutron star binaries. Although a magnetic mountain is 
susceptible to transient, toroidal, ideal-MHD instabilities, 
these are not disruptive. The saturation state still confines 
the accreted matter to the magnetic pole, efficiently screen- 
ing the dipole moment in the long term. 

This article is concerned with the fate of a magnetic 
mountain when a nonzero electrical resistivity switches on. 
We extend the ideal-MHD code zeus-mp to add a resistive 
term to Ohm's law and perform three-dimensional simula- 
tions for different values of the resistivity. In the axisym- 
metric case, we find that global MHD oscillations compress 
the magnetic field, accelerating plasma slippage across flux 
surfaces. As a consequence, the mountain relaxes on a time- 
scale ti which is shorter than the diffusion time-scale td but 
comparable to the accretion time-scale. In the nonaxisym- 
metric case, Ohmic diffusion additionally tends to restore 
axisymmetry. We do not find any evidence of transient re- 
sistive instabilities, like the resistive ballooning mode, on the 
intermediate tearing mode time-scale (tatd) 1//2 . The moun- 
tain persists over ~ 10 5 — 10 8 years, comparable to the du- 
ration of the accretion phase in a low-mass X-ray binary 
(LMXB). 

Astrophysically, the key result of the paper can 
be stated as follows: magnetically confined mountains in 
LMXBs are stable (in ideal and nonideal MHD) over the ac- 
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cretion time-scale and relax over the typical life-time o f radio 
millis econd pulsars (~ 10 9 yr) after accretion stops. Ijonesl 
(|2004l ) argued that the electrical conductivity in the solid 
crust is significantly lower than that for a homogenous bcc 
lattice and temperature- independent, with r\ — 1CP 24 s. For 
this value, we expect a stationary state at M a ~ 10~ Mq 
where the diffusive mass flux escaping the polar cap is ex- 
actly replenished by accretion. To study the structure of 
such a state, and confirm its existence, we must extend the 
growing simulations in Fig. 1121 a key topic for future work. 

One shortcoming of the calculations is the neglect of 
rotation. Ac creting millisec o nd pulsars spin up as fast as 
fi ~ 620 Hz (|Gallowavll200l ). ISpitkovskv et all (|2002h found 
that surface thermonuclear burning is unaffected by rotation 
in its early stages, but the thermonuclear flame spreads more 
slowly as time passes. iBhattacharvva fe Strohmaverl (|2007T ) 
applied this idea to qualitatively reproduce the light-curves 
from 4U 1636-536 and SAX J1808.8-3658. The Coriolis 
force also modifies the continuous part of the ideal-MHD 
spectrum for ax isymmetric configurations with a uniform 
angula r velocity l|Hellsten fe Spiesll 19791 ; IVigelius fe Melatosl 
l2008al ). especially for short- wavelength modes. However, it 
does not affect the equilibrium configuration. A simple es- 
timate shows that one requires a transversal speed of v ~ 
8 x 10 8 cm s" 1 to attain a Coriolis forc e which is compara - 
ble to c 2 VP. The star may also precess (|Chung et al.ll2008l ) . 
complicating the treatment of rotational effects. 

Magnetic mountains in LMXBs are promising 
sources of gravitation al waves l|Melatos fe Payne! 120051 ; 
iPavne fc Melatosl 120061 ). For M a > 1O" 5 M , there is a fair 
prospect of detection with next generation interferometric 
detectors like the Laser Interferometric Gravitational Wave 
Observatory (LIGO) (VM08). Clearly, the strength of the 
signal depends critically on the long-term stability of the 
mountain and the rate at which it relaxes res istively. In a 
companion paper (|Vigelius fc Melatosl l2008bh . we predict 
the signal-to-noise ratio attainable by LIGO when the 
resistive results of this paper are included. We also show 
that the electrical resistivity can be constra ined by existing 
LIGO data, by invoking the lBildstenl (|l998l ) torque-balance 
limit for LMXBs and th e Blandford spin-d own limit for 
radio millisecond pulsars (|Abbott et al.ll2007h . 
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APPENDIX A: IMPLEMENTING RESISTIVITY 
IN ZEUS-MP 

Al Advection step 

Resistive MHD comprises a set of seven coupled, nonlinear 
partial differential equations for the magnetic field B, the 
bulk velocity v, the plasm a density p, and the pressure p 
jGoedbloed fc Poedts|[2004l ) : the equation of mass conserva- 
tion, 

| + v.(H = o, 

the momentum equation, 



(Al) 



/ dv 



p + v ■ Vvj + Vp - (V x B) x B + pVifi = 0, (A2) 
and the induction equation, 

^-Vx(»xB-ijVxB)=0, (A3) 

where r\ denotes the resistivity and ip is the gravitational 
field. The system is closed by the supplementary condition 
V-S = and an isothermal equation of state p — c^p, where 
c s represents the isot hermal sound speed . 

As explained by lHaves et al.l (|2006h . ZEUS-MP employs 
an operator split algorithm based on the method of finite 
differences on a staggered grid. The advection step is done 
in two stages. Firstly, a source step solves 

P^ = -V P -V Q pV^-V{B 2 /2p Q ), (A4) 

where an artificial viscous pressu re tensor Q is in- 
cluded (|von Neumann fc Richtmverl Il950l ). Secondly, to 
treat transversal MHD waves properly, one must solve the 
magnetic tension force along with the induction equation in 
a single step using the metho d of characteristics and con- 
strained transport (MOCCT) (|Hawlev fc Stonelll995l ). The 
MOCCT step advances B by computing the line integral 
of the electromotive force (EMF) e = v x B around a cell 
boundary S: 



dt 



B -dS 



edl. 



(A5) 



Second-order accuracy in time is achieved by employing 
time-centered values for e. The extrapolation in time is 
done using the characteristic equation for transverse Alfven 
waves. It can be shown that MOCCT ensures V • B = 
to machine accuracy at all times provided the initial field is 
solenoidal. The extrapolated B is then used to work out the 
transverse magnetic forces and accelerate the fluid accord- 
ingly: 



dv 



dv 
dt 



sourccstcp 



- / i ( 7 1 (B-V)B. 



(A6) 



Finally, the fluid density and momentum are advected via 



df 
and 



dt 



pdV = 



pv dV 



pv-dS 



pvv ■ dS. 



(A7) 



(A8) 



OY 



A visual comparison of ()A3|) and (|A5I) suggests a natural 
way to incorporate the resistive term (|Stonelll999l ): we use 
the updated B to work out the current density j = V x B 
and apply equation (|A5|I again, replacing e by —r\j this 
time. The staggered grid allows for central differencing and 
thereby guarantees spatial second-order accuracy. This re- 
sistive algorithm has been used in conjunction with ZEUS-3 D 
to study protostellar jet formation fFc ndt fc Cemeliio!l2002l ). 

A von Neumann analysis in terms of eigenmodes yield s 
a stability criterion for parabolic PDEs (|Press et alj|l986h . 
At ^ 2A 2 ?7~ 1 , which depends quadratically on the minimal 
grid cell size A. We find empirically that our implementation 
requires At ^ 10 _2 A 2 r; _1 . If 77 is high, this constraint dom- 
inates the ideal-MHD timestep and drastically increases the 
run time. We therefore make u se of a superstep algor ithm, 
similar to the one described bv lAlexiades et al] (|l996l ). We 
compute the ideal-MHD timestep AfcFL according to the 
usual Courant-Friedrichs-Levy (CFL) condition, as well as 
the resistive timestep Assistive- After updating B by the 
MOCCT procedure, we apply the resistive algorithm in a 
cycle of N steps, such that A£ C fl = AAT, with AT < 
Afresistive • This approach was implem ented successfully in a 
resistive module for the pluto code |Mignone et al.|[2007h . 

Special care must be taken when incorporating the 
boundary conditions, zeus-mp adds two and three ghost 
cells at the inner and outer boundaries, respectively, where 
it either sets e according to the boundary conditions for B 
and v or communicates it at a processor boundary (inside 
the integration volume) via the message passing interface 
(MPI). The processor boundaries are set by the MPI topol- 
ogy, i.e. the division of the computation grid among the dif- 
ferent processors. The staggered grid requires e at the inner 
boundary, so we need to add another layer of ghost cells for 
B and v there to provide e in all ghost cells. In order to 
minimize the alterations to the code, we prefer to compute 
(or communicate to the neighbouring processor) the whole 
layer of ghost cells for B at the beginning of every super- 
step. We employ t he MPI communication flow described in 
lHaves et all (|2006ft to minimize inter-processor traffic. 

A 2 Test case 

We test our code extensions by simulating a purely diffusive 
problem. We set p = 10 9 po (see section I2.1|l and v — 
to suppress any fluid motions. Equations <| A If) (IA3|) then 
reduce to a single diffusion equation 

OB 
~dt 

which can be solved easily in Cartesian coordinates (x, y, z): 
— < A1 °) 



-V x (77V x B), 



(A9) 



B(r,t) = 
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Figure Al. Snapshots of the three-dimensional diffusion problem in spherical coordinates at t = to (stars) and f = to + 2t 1 { (crosses), 
where = 0.22to, along with the analytic solution (I A 1 H — (I A14 1 (solid curves). Shown are B r , Bg, and (left, middle, and right panels 
respectively) for three cross-sections: r = 5/iq, 8 = 1.18 rad, and cf> = 3.14 rad. The agreement is excellent (2.4 per cent). 



r -(j/ 2 +z 2 )- . -(ai 2 + z 2 )- . -faj^+y^l » i 

x|e K " 'e x +c K 'e y +e 1 " e. z \. 

A coordinate transformation yields the result in spherical 
coordinates: 



B r (r,t) 



t 

x ( e 



-2 ,™2 



I) 



e 4 <n cosfl + sin6» (All 

™ 2 ^j!, _a .,„2 



4?7t 



cos <p + e" 



tin^ 6 sia~ 4> 



477* 



sm < 



»(*•>*) 



e ^ 



(A12) 



r 2 cob 2 sin 2 fl 



- 2 .;„ 2 o.;„ 2 . 



x e 4 i' cos q> + e 4 'j* sm ( 



and 



-e 4t <* sm t 



e 4 i* 



(A13) 



(A14) 



r 2 sin 2 8 sin 2 <j> 
x I e 4 '!* COS ( 



t* 2 cob 2 <t> Bin 2 8 

e 4t )' sin < 



Test runs were performed in Cartesian (three dimensions) 
and spherical polar coordinates (two and three dimensions) 
employing time-dependent boundary conditions with Eqs. 
(|A11|) and (|A14|) . Fig. |A1 shows two snaphots of the three- 
dimensional run in spherical polar coordinates at t — to 
(stars) and t — to + 2rd (crosses), with td = 0.22ro, along 
with the analytic solution (|A11|| - (|A14|) . They are in ex- 
cellent agreement, with relative error < 2.4 per cent at 
t = To + 2r d . 



APPENDIX B: RESISTIVE PARKER 
INSTABILITY 

In this section, we derive an analytic dispersion relation for 
the linear, resistive , MHD modes of a pla ne-parallel, gravi- 
tating plasma slab (|Singh fc Tandonll 19691 ). The ideal-MHD 
cou nterpart of this problem is known as the Parker instabil- 
ity (|Parkedll967l : lMouschoviaslll974l ). 

Let us assume a uniform gravitational acceleration g di- 
rected parallel to the z-axis, and a unidirectional magnetic 
field parallel to the y-axis. We can then write down the mag- 
netostatic equilibrium. The density and magnetic field are 
given by 



po exp 



-'/- 



u 2 (l + a) 



B(z) — Bo exp 



-gz 



2v?{l +a) 



(Bl) 



(B2) 



under the additional assumption that the magnetic pressure 
is proportional to the gas pressure everywhere, viz. 

The equation of state is p — u 2 p. 

Next, we write down the linearized equations of mass 
conservation, 



at 

force balance, 

dv<» 
P- 



dt 



- Vp W _ J_ V \2B ■ B w ] 

2p L J 



(B4) 



(B5) 
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_l_ _L^J . VB*- 1 ' H — B^ 1 ' ■ S7B — p' 1 'ge,(B6) over the time required to diffuse across one hydrostatic scaie 

Mo Mo height. 

and induction, 

= l^B 11 ' + V x \ V W x Bl . (B7) 
at 4:71(7 1 J 

In |]H])-(|Bg and below, p (1) , p (1) , B (1) , and i> (1) denote 
the perturbations of the pressure, density, magnetic field, 
and velocity respectively, a denotes the conductivity. 

Ignoring interchange modes (k x — 0), we assume the 
perturbed quantities have the form oc e ^ k vv- ut ) _ Further- 
more, we only consider perturbations in the y-z plane, 
i.e. v x 1] = and B w = [0, dA x l) /dz, -dA^/dy], where 
A = A x x is the vector potential. Eq. ()B4|I yields 



- ilu P w + ikypvW + p-^L - l v Wp = o, (B8) 
oz L 

with L = m 2 (1 + a)/g. Similarly, the components of ()B5[) 
reduce to 

- 1UJ pvM = -m%pM+BU^, (B9) 

Ltty 



and 



(i) 2dp {1) dA ( ^au 2 p 
iujpv y z = -u — — - (BIO) 

OZ OZ LBy 



By 

fJ-0 



d 2 ^ d 2 A x 1} 



P (1) g- (Bii) 



dy 2 

Finally, the induction equation (|B7|I yields 

In order to solve (|B8[) - (|B12[) analytically, we make the 
short- wavelength approximation d/dz <C k y . Eliminating 
v* and p' 1 ', we find the dispersion relation 

+ 2^kl _ 2au 2 ky] (lu 2 - u 2 kl) = J^L. (B13) 
For a -»• do, (|BT3| can be solved to obtain 



2u 2 <J = (1 — 2a)u 2 k 2 ± 



(l-2a)Vfc*- 4 ^ Q 



(B14) 



1 + a 

The modes are stable when the discriminant is positive: 

!f!*» > 4q 1 fB 15) 
ff 2 > l + a(l-2a)2' (B15) 

For large but finite a, lu x is perturbed slightly, with lu — 
cooo + lu' and \lu'\ <IC Woo- Solving for lu' , we obtain two 
branches, the first damped, 

i 2 2 

^damped = 5 = — ~ (B16) 

07TCT ZTD 

and the second growing, 

, 2 







167TO- 1 + a (1 - 2a) 2 w 4 

1 Q g 2 i 

2(1 - 2a) 2 l+ak 2 u 2 To' 



(B17) 
(B18) 



The damped mode has a decay time roughly equal to the dif- 
fusion time td , whereas the growing mode amplifies quickly, 



